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Continuum theory of memory effect in crack patterns of drying pastes 
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Department of Applied Mathematics and Physics, Tottori University, JP-680-8552, Japan 

(Dated: April 7, 2008) 

A possible clarification of memory effect observed in crack patterns of drying pastes [A. Nakahara 
and Y. Matsuo, J. Phys. Soc. Japan 74, 1362 (2005)] is presented in terms of a macroscopic elasto- 
plastic model of isotropic pastes. We study flows driven by steady gravitational force instead of 
external oscillation. The model predicts creation of residual tension in favor of cracks perpendicular 
to the flow direction, thus causing the same type of memory effect as that reported by Nakahara 
and Matsuo for oscillated CaCOa pastes. 



I. INTRODUCTION 

As the plastic behavior of soft dassy materials has 
been attracting increasing interest [l|, it was reported 
by Nakahara and Matsuo [2, S |j] that a drying paste 
exhibits a memory effect. They observed a drying pro- 
cess for a paste containing calcium carbonate (CaCOs) 
and water in a shallow container in order to study the 
resulting crack pattern. The crack pattern was typically 
found to be isotropic, but they discovered a way to intro- 
duce anisotropy into the paste before the drying process 
commences: by applying a horizontal oscillation to the 
container immediately after the paste is poured into it, 
a memory of the oscillation is imprinted into the paste, 
which determines how it should break in the future. 

Through systematic experiments, Nakahara and Mat- 
suo also found that plasticity is essential to the memory 
effect in CaCOs pastes. No memory effect is observed 
if the strength of the applied oscillation is below the 
threshold value corresponding to the plastic yield stress 
of the paste. Just above the threshold value, the paste 
remembers the oscillation that caused the plastic flow, 
developing cracks perpendicular to the direction of the 
oscillation. If the oscillation is too strong or the paste 
contains too much water, waves and global flows are in- 
duced, eliminating the memory effect. A different kind of 
paste (mixture of magnesium carbonate hydroxide with 
water) l5| exhibits not only a memory effect similar to 
that of CaCOs that occurs just above the threshold of 
plastic flows and causes cracks perpendicular to the ex- 
ternal oscillation, but also a different type of memory 
effect in its water-rich condition where the cracks are par- 
allel to the direction of the global laminar flow caused by 
the oscillation. Too strong an oscillation and too much 
water also destroy the memory in this paste, with the 
emergence of chaotic, turbulence-like flows U^ charac- 
terized by fluid motion in every direction. 

Here we focus our attention on the former type of mem- 
ory effect that causes cracks perpendicular to the external 
forcing, which we refer to as the Type-I Nakahara effect. 
The latter type, which could be called the Type-II Naka- 
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hara effect, will be discussed only briefly. 

Although it is certain that the memory effect in CaCOs 
pastes originates from plastic flow, it is unclear which as- 
pects of the plastic flow are essential. More specifically, 
because the role played by the unsteadiness of the flow 
is not fully understood, it is unknown whether a slope 
flow, in which the external forcing is steady, can cause 
a memory effect. To answer this question, coworkers of 
the present author have started an experimental study on 
the slope flow of CaCOs paste. Paste supplied through a 
funnel is driven downstream by gravitational force, and 
when the supply is stopped, the paste "freezes" at some 
finite thickness due to the finite yield stress. Preliminary 
results suggest the presence of a memory effect (Type-I 
Nakahara effect), where the cracks are perpendicular to 
the direction of the flow, i.e. the direction of the exter- 
nal forcing. Details of the experiment will be reported 
elsewhere \&\. 

As a first step in the theoretical investigation into the 
slope flow of CaCOa paste, we study the dynamics of an 
elastoplastic liquid layer with constant thickness falling 
down an inclined wall. First, we construct a continuum 
model equation that meets several requirements, so that 
it can be a good description of the CaCOs paste. Next, 
we apply this model equation to the two-dimensional 
slope flow with constant layer thickness. We find that the 
flow develops tension in the streamwise direction, which 
remains in the paste. Since the residual tension implies 
that the dried paste will be more fragile in the pertinent 
direction, this result presents a possible clarification of 
the Type-I Nakahara effect. 



II. REQUIREMENTS FOR THE MODEL 

The strategy in this paper includes the construction 
of a set of model equations acceptable as a continuum 
description of CaCOa paste. A useful precedent for this 
model construction can be found in the continuum me- 
chanics of gases and simple liquids [T], [g], in which the 
Navier-Stokes equation is deduced from several macro- 
scopic requirements, such as homogeneity, isotropy, and 
the postulation that the deviatoric stress tensor is a lin- 
ear function of the rate-of-strain tensor (without time 
lag). Following this precedent, let us list the analogous 



requirements for paste flows. 

We assume that the dynamics of the paste is isotropic, 
in the sense that the paste has no preferred direction ex- 
cept for the principal axes of the stress tensor. This is 
plausible for CaCOa, which consists of spherical parti- 
cles 0, 01- On the other hand, magnesium carbonate hy- 
droxide is not expected to exhibit isotropy in this sense, 
as its particles are disk-like |^5] and therefore can exhibit 
anisotropy similar to that of liquid crystals. 

The stresses in the pastes under present consideration 
are primarily sustained by the interparticulate bond net- 
work. There should be also a contribution from the vis- 
cosity of the solvent (water), but we assume that this 
contribution is much smaller than that of the interpartic- 
ulate bonds (in other words, we consider only very thick 
colloids). Unlike the chemical bonds, the interparticu- 
late bonds in flowing pastes are usually so breakable that 
they are constantly destroyed and reconstructed. The 
stress is therefore expected to be governed by a Maxwell- 
type equation [3, Uo, [ill, 113, [3 whose relaxation time 
represents the lifetime of the bond. 

We postulate that the relaxation time, denoted by r, 
is a scalar: the collapse of the force network involves 
bond breakage in all directions. Since the paste is plas- 
tic, the relaxation time r must be variable. An infinitely 
large r represents solid- like behavior, while a finite r de- 
notes fluidity. The transition between these two behav- 
iors with a certain threshold gives a formulation of plas- 
ticity. Isotropy dictates not only that r itself is a scalar, 
but also that r should be a function of some scalar quan- 
tity. With the von Mises criterion [lj| and its energetic 
interpretation [1^ in mind, we assume that the relax- 
ation time r is a function of strain energy. Introducing e 
to denote the nondimcnsionalized strain energy (defined 
later), this assumption is formulated as 



r = r(e) 



+00 (e < threshold) 

'''0 = 7?p/5' (£ ^ threshold) 



(1) 



where r\^ is a constant with the dimension of viscosity, 
and S is the shear modulus. 

We describe the relaxation of the bond network in 
terms of the Lagrangian (material) variable ^, rather 
than the Eulerian variable r. The main reason for this 
choice is the adequacy of the Lagrangian description for 
tracing so-called frozen quantities. With the relevant 
physical quantity provisionally symbolized as Q (prob- 
ably representative of the density of the bond network) , 
the equation of relaxation is expected to have the form 



1 



^|K(tt) 



Q.ii.t)- 



(2) 



In the limit of an infinitely long relaxation time (t -^ 
+oo), Eq. ^ reduces itself to 

dtQ{i,t) = ^ (3) 

which manifests directly that Q is "frozen" in the mate- 
rial. In the Eulerian description, the same assertion as 



Eq. ([3]) would have a more complicated form, 

(9t -I- V • V) g(r, t) + • ■ • = (Eulerian), (4) 

where "• • • " stands for various convective terms required 
according to the tensorial character of Q . Since we pre- 
fer the clarity of Eq. ([3]) to the obscurity of Eq. (|4]) , the 
Lagrangian description is adopted during the construc- 
tion of the model (the result could be reformulated in 
the Eulerian description after it is developed, of course) . 
Generally, the mathematical formulation of elasticity is 
related to the deformation of the fiuid (or material) ele- 
ments. In the steady and quasi-steady motions of pastes, 
deformation (as opposed to rate of deformation) can in- 
crease unlimitedly as the time elapses. This requires our 
model to be free from the restrictive assumption of a 
small deformation, motivating the inclusion of the geo- 
metrical nonlinearity to the full extent. Besides, the solid 
behavior of the paste for a small deformation should have 
an isotropic Hookian limit (with shear modulus S), be- 
cause the paste is isotropic. For the same reason, the 
fluid behavior is expected to have a Navier-Stokes limit 
for small r or small shear rate (which is realized for water- 
rich pastes with a vanishingly small yield stress) . In both 
behaviors, we regard the paste as incompressible, as far 
as flow processes are concerned, neglecting the slow ef- 
fects of drainage and evaporation. Finally, the model 
equation must have "relabeling symmetry" [16], i.e. the 
system of equations must remain formally unchanged in 
regard to the change in the Lagrange variables. In what 
follows, while making some additional assumptions, we 
will construct a system of model equations that satisfles 
all of these requirements. 



III. MODEL 

In this section, we construct a continuum paste model 
for a generic ri(j-dimensional geometry. The model equa- 
tions will be summarized at the end of ijIIIBl Subse- 
quently, in i jIVI and i fVl this model will be analyzed under 
a specific setup describing a two-dimensional slope fiow 
with constant layer thickness. Readers who are more in- 
terested in the analysis than the model construction may, 
after checking Figs. [T] and [H skip to Eqs. (M)) at the end 
of gWl 



A. Kinematics 

First, we review the Lagrangian description of kine- 
matics. The configuration of an nd-dimensional contin- 
uum is represented by a mapping from Lagrangian vari- 
able ^ (also known as "label" or "material variable" [Ifl] ) 
to the position vector r. For rid = 3, we write 



C=(e,^,C)^r = r(|,t) 






(5) 



where 



denotes the representation in terms of Carte- 



sian components. For n^ = 2 we will omit 77 and y, as- 
suming that all the motion occurs in the (x, z)-plane. 
The time-derivative of r = r{^,t) gives the velocity, 



^r = dtvitt). 



(6) 



In Eq. ^ and in what follows, dt stands for the 
time-derivative in the Lagrangian description (Lagrange 
derivative, which is usually denoted by D/Dt in Eulerian 
description). Using {(9^r, 9^r, 9^ r} (where di = d/d£,^) 
as the set of local bases, we can represent the velocity as 



V — v^diV. 



(7) 



In Eq. ([7]) and in what follows, summation over i E 
{Cj Vj C} is understood according to Einstein's contraction 
rule. The coefficients (v^) in Eq. ([7]) are referred to as the 
contravariant components of v (see Eqs. (|A1[) and (|A6|) 
in Appendix El . The acceleration is 9tv = dt (ti'^^r); we 
emphasize again that dt denotes the Lagrange derivative. 
The square of the Euclidean distance between two 
neighboring "particles," labeled by ^ and ^ + d^, is 

ds^ = \{d,r) def = 9^Jded^', 9^J = (5,r) • (5,r), (8) 

which introduces the metric tensor denoted by {gij) or g. 
In this paper we refer to g as the "Euclidean" metric ten- 
sor, which does not mean that Qij is equal to Kronecker's 
delta but means that the Euclidean metric of the r-space 
is imported into the ^-space by Eq. ([5]). 

In general, it is totally unnecessary to choose ^ to be 
some "initial" position of the element, except for some 
particular situations in which the initial state has a spe- 
cial significance. One of these special cases is that of 
purely elastic bodies initially set in a stress-free and unde- 
formed state, called a "natural state" [13| . It is meaning- 
ful in this case to choose the "natural state" position vec- 
tor as ^ so that gij defined by Eq. ([5]) is essentially iden- 
tical to the Cauchy-Green deformation tensor l9] whose 
difference from 5ij is responsible for the elastic restoring 
force. This is a rather special case, however. More gen- 
erally, ^ has nothing to do with the initial state, and the 
natural metric tensor ^ is used as a reference to define 
the elastic deformation, instead of assuming the global 
existence of the stress-free natural state. The (locally) 
undeformed state is formulated as g = g'' , and the differ- 
ence between g and g^ is responsible for the stress. More 
details about g'' will be discussed later. 

The incompressibility condition is expressed as 



dt det g == 0, 



(9) 



because the mass of a fluid element is pVdetgd"'^^ which 
should remain unchanged, and the density p also remains 
unchanged during the motion. For simplicity, we assume 
that p is a global constant. Then, without loss of gener- 
ality, we can replace Eq. ([9]) by 



detg= 1. 



(10) 



B. Equation of motion and constitutive relation 

Now we detail the dynamics. With the stress field 
denoted by P = P^^{dir) (g) (djv) and the external body 
force by F = F^diV, the momentum equation is written 
as 



p dt {v'd^r) 



d_ 
dr 



'{P''{d^r)<»{djr))+F'd^v 



or, in contravariant component representation, as (ij 



P{d, 



'tV + w-'VjU* 



-V,P^- 



F' 



ill) 



The left-hand side is the contravariant component of 
the acceleration vector 9tv multiplied by the density p, 
and Vj denotes the covariant derivative (these math- 
ematical concepts are clarified in Appendix [X] to the 
degree sufficient for the present work; for a more pro- 
found understanding of the mathematical background, 
see Refs. [Ill, [l8|)- While F is regarded as given, P must 
be determined by a suitable constitutive relation. 

From the discussion in the previous section, we expect 
that P obeys a viscoelastic equation of Maxwell type. 
The Maxwell model is often illustrated as a spring and 
dashpot connected in series Q, for which the relation 
between the tension T and the total length x is given by 



r = 



K {xs - xi] 



dXY) 



X — xs + xy) (12) 



where k is the spring constant, p is the resistance, a;g and 
xd are the length of the spring part and the dashpot part, 
respectively, and Xg denotes the natural length of the 
spring part. It is customary to eliminate the "internal" 
variables (xg, xd and Xg) from Eq. (|12p . which yields 



M 



-'s'- 



dx 
di' 



(13) 



A timescale p/ k in regard to stress relaxation is recog- 
nized in Eq. (fH)) . 

Now it is necessary to elaborate the Maxwell model in 
two respects: it needs to include plasticity and it also 
needs to describe n<j-dimensional continuum mechanics. 
In regard to the first point, most of the existing studies 
are based on an elasto-plastic decomposition, which is 
a direct extension of Eq. (H^). However, this approach 
has a disadvantage in that the incautious use of internal 
variables can lead to a difficulty, in particular for a finite 
deformation JJ[, _20] . Here we adopt a different approach 
that is closer to Eq. P^ . thereby avoiding a direct ref- 
erence to the internal variable xy>- 

The essential idea is to attribute the relaxation to the 
natural length x^ , which is related to the tension T as if 
the model is totally elastic: 



(14) 



T = k{x — x^). 

The natural length x^ can be expressed as x'' = a; 

in terms of internal variables, but this relation is not to 



^D 



be used explicitly; we note only that a;'' is time-dependent 
while Xg is not. By substituting Eq. pi)) into Eq. P^ . 
we find an equation that describes the relaxation of the 
natural length a;'' : 



— — = — (x — x^ 
di n ^ 

or, by introducing r = /i/^, as 

l + T—-\x^=x 
at 



(15a) 



(15b) 



in the form of relaxation toward x^ — x. Eqs. (|14p and 



P^ provide us with a prototype of the plastic model. 

Let us find nd-dimensional continuum equations corre- 
sponding to the prototypical equations (I14p and (I15p . As 
a candidate, we adopt an elastic constitutive equation 



together with an inelastic equation 



9, 



't9^ 



V _ 



-'^s^ 



v*g- 



(16) 



(17) 



where {g^^ ) denotes the inverse of the component matrix 
of the "Euchdean" metric tensor (g^), and (gl^) is that 
of the natural metric tensor, such that 



9ij9^'' = 5, 



u5t, - <^^ ■ 



The natural metric tensor g^ represents the square of 
the "natural distance" between two neighboring points 
labeled by ^ and ^ + d^. 



(di 



gh'^Cde 



(18) 



in the sense that the difference between ds^ and fds'') 
accounts for the restoring force according to Eq. ([16]) . In 
the special case of purely elastic bodies initially set in a 
stress- free "natural state" (at t = to), ds^ is the distance 
in this initial configuration and g'' is the corresponding 
metric: 



9i 



(purely elastic case) 



In general, however, g'' differs from the initial value of 
g. This is inevitable due to Eq. (|17p . which prescribes 
that the natural metric g'' is subject to relaxation. To 
make Eq. (fT7|) more easily recognizable as a relaxation 
equation, we rewrite it in the manner of Eqs. ([2|) and 
(|15bp as 



{l + Tdt)g;'=Kg^^, T 



K 






(19) 



this equation provides that {gV ) should evolve toward an 

isotropic tensor {Kg^^). Plasticity is incorporated via r 
according to Eq. JT]) . The idea of using a natural metric 



to reformulate the Maxwell model has been known among 
several researchers of rheology (including the authors of 
Refs. [lil[l3), but the present author could not identify 
any publications in which the notion of the natural metric 
and its relaxation is formulated explicitly. 

The nd-dimensional elastic equation (I16p . correspond- 
ing to the one-dimensional Hookian equation (fT4|) . orig- 
inates from consideration of elastic strain energy. Since 
gijd^'^d^^ is a positive definite quadratic form, there ex- 
ists a set of Euclidean vectors {P5,Pr;,Pc} such that 
9ij — Pi ■ Pj (this is proved essentially in the same way 
as the polar decomposition theorem [3, [l3|)- Then, by 
defining 



dr^ = p, de 



(20) 



we have (ds^) = dr^ • dv^ = |dr^| . Note that Eq. ^ 
does not claim that dr'' is a differential of "r''": such in- 
tegrability is not guaranteed. However, it is legitimate to 
interpret dr^ as a natural configuration of each small ele- 
ment. Since p^'s must be linearly independent due to the 
positivity of detg'', Eq. (PO)) can be inverted, which we 
denote as d^* — p* • dr'' . From this and the "Euclidean" 
metric (jS]) , we have a relation between the Euclidean dis- 
tance ds and the natural configuration dr'' , 



ds2 = (g.^. p*^ ® Pi) : (dr^ ® dr^ 



(21) 



Let us denote the eigenvalues of this quadratic form by 
{A^} so that ds^ — \\ \dv^\ along the a-th principal 
axis. The geometrical meaning of Aq is clear: it repre- 
sents the elongation factor of the line element. Isotropy 
requires that the elastic energy (denoted by E) should 
consist of a symmetric combination of these eigenvalues. 
The simplest form with a correct Hookian limit is 



E=\s{\l + \l + \l-i) 



(22) 



for rid = 3. E g. (1^^ is known as neo-Hookian constitu- 
tive equation 17] . With the aid of the incompressibility 
condition, which implies A1A2A3 — 1, Eq. (j22p reduces to 



E = S{ei + 62 + ef) for small deformations (Aa = 1 + ec, 
and JCqI <C 1). By using the definition of p* and in- 
troducing e = J2a i^a ~ 1) 1 foi' finite deformations, the 
elastic energy E is expressed in terms of the inverse nat- 
ural metric tensor: 



E^\se. 



9ij9l^ -rid- 



(23) 



By calculating the variation of the elastic energy E 
in regard to r through the metric tensor g under the 
constraint of incompressibility condition pop . we find 
that the contravariant components of the stress tensor 
are given by Eq. p^ . Details of this calculation are 
shown in Appendix [Bl Note that the tensor (5*^) in the 
first term of the right-hand side of Eq. P^ stands for 
the Euclidean unit tensor: 



g'^ {d,r) (g, (djr) = 1. 



(24) 



Thus we find that the term pt/'^ stands for an isotropic 
stress. The scalar p is related to the hydrostatic pres- 
sure arising as a constraint force (Lagrange multiplier) 
for incompressibility. It is convenient to define 



S [g;^d,r) ® id,r) 



1 



(25) 



and call it the "elastic stress tensor" so that the stress 
tensor P = P^^ {div) ® {djv) is given by 



pt 



(26) 



It is easy to confirm that ct vanishes when g^, — gij . 

We emphasize that {gl-' ) in Eq. ((25|) , which determines 
the elastic stress tensor tr, is the inverse of the natural 
metric tensor. This must be the case so that a should 
remain invariant under the relabeling of the Lagrange 
variables. This is also acceptable if we remember that 
springs with different lengths but the same local proper- 
ties obey a constitutive relation analogous to Eq. (|^5|) . 



r-^o(^-l) 



where Sg is the normalized spring constant, and notice 
that a is an intensive variable as well as the tension T 
and therefore must be expressed as such. 

Now we discuss the inelastic part of our model de- 
scribed by Eq. (fl7|) . This equation states the relaxation 
of the inverse natural metric tensor, formulated accord- 
ing to the following discussion on interparticulate bonds. 
The natural metric represents the energetically optimal 
configuration of the particles determined by the bond 
network. The network strength, or the bond density, is 
represented by the inverse natural metric g\^ (not by the 
natural metric g' itself). In flowing pastes, however, this 
bond network is ephemeral. We suppose that the network 
is destroyed at some rate and reconstructed isotropically. 
With the destruction rate denoted by ly and the recon- 
struction rate by ;/,, the temporal change of the bond 
density is given by —i^gl'' + J^* <?*-', leading to Eq. (fT7|) . 

The ratio K = v^:/v is determined by postulating the 
incompressibility of g^ , 



detg^ = 1. 



(27) 



Differentiating Eq. (I27p with regard to t and then sub- 
stituting Eq. (J17p into it, we find that n^v = g^^g^-'i^*, 
which implies 



K = 



rid 



rid 
9%9'' 



(28) 



According to Eq. ([T]) in the previous section, r is sup- 
posed to be a function of the elastic strain energy, so 
that T — T(e) with e given by Eq. (P5)) . The simplest 
form consistent with Eq. ^ is 




(o^/sf 



FIG. 1: Tlie inverse of tlie relaxation time, r~^ — ^{e), 
defined by Eq. (|29[) . According to neo-Hookian constitu- 
tive equation of nd-dimensional elastic bodies, e is given by 
Eq. ^. 



where ay is the yield stress (we will see later that the 
energy e for shear stress cr is calculated to be a'^/S^). 
It is physically more realistic and mathematically less 
problematic to suppose that r is a continuous function 
of e. Here we assume 



He) 



1. 



max 0, 1 



o-y/S 



(29) 



which is Lipschitz-continuous in spite of weak singularity 
at the yield point (Fig. [T]). Eq. ([29]) is chosen in such a 
way that it agrees with Bingham plasticity [l^ [2l|, [22| for 
simple shear flow with shear rate 7, where the shear stress 
a is estimated to be cr ~ St'^. Admitting e ~ a^/S^, 
from Eq. ((29l) we find 



— = i'{e) Ve = 







'{S^-ay) 



[e < 4/5^) 
{e > 4/5^) 



which is Bingham plasticity. 

Let us summarize our model. The governing system 
of equations consists of Eqs. HI]), ^, HH), and ((29|) . 
supplemented with the kinematic relations ([5]), ([7]) and 
([8]), as well as incompressibility conditions pO)) and (|27p . 
Eq. (I29p requires the evaluation of e by Eq. ((23|) , which 
is actually not independent of Eq. (flB]) . but should be 
included in the model for convenience. The independent 
variables are ^ and t (Lagrangian description), and the 
essential dependent variables are r and g''. The velocity 
and the Euclidean metric tensor are derived from the 
differentials of r — r{^,t). Due to the incompressibility 
condition, there arise two additional scalar fields, namely 
p and K; the latter is determined by Eq. ((28|) . 



C. Navier-Stokes limit 

There remains the task to confirm that the whole sys- 
tem of model equations reduces to the nd-dimensional 
incompressible Navier-Stokes equation if t is set to be a 
small constant such that r <C ||Vv||~^. By expanding g- 
(as well as gl-' ) and K in power series of r, from Eqs. (fT9|l 
and (|27p we find 



91^-9'' 



rdtg^^+0{T^), 



K = 1 + 0(t2). 



(30) 



z=H 



z=0 




FIG. 2: Schematic view of the system and coordinates. A 
uniform fluid layer with thickness H is assumed. The la- 
bel variable ( coincides with the depthwise Eulerian coordi- 
nate z in the present setup, and the velocity is Usx where 
U = dtX{(,t). The gravitational acceleration vector is G = 
{Gsm9)ex — (G cos 6)62:, whose x-component is Gx = Gsin^. 



The time-derivative term dtg^-' on the right-hand side of 
Eq. ([50]) is calculated as 



dtg'' = -g'^Q^'d, 



hgki 



(31a) 



and 



^tg^J = dt ((ftr) • {djY)) 

= (9,v) • {d,v) + (a,r) • (a,v) 



(31b) 



where {vi) denotes the covariant components of the ve- 
locity vector V, and V^Wj denotes the covariant derivative 
oivj defined by di{vjVS,^) = {VrVj)V£,K Using Eqs. jST]) 
to evaluate g!"' in Eq. ([30)) . from Eq. (|25|) we obtain 

a^S (g;={d,v)®{d,v)-t) 

= S {g^ - g'') (a.r) ® {d,v) 
= -St {dtg'') {d,v) ® (9,r) 
= Sr [g'^g^'dtgu) {^^v) ® {d,v) 
= ST{dtgki){Ve)®{^S!) 

= St (V(8)v-H'(V(g)v)) , 

and identify it with the Newtonian-Stokesian relation 

(7 = 277*symgradv (77* — St) (32) 

where symgradv denotes the symmetric part of V v. 
The equation of motion then reduces to the Navier-Stokes 
equation, which was to be demonstrated. 



IV. SIMPLIFICATION FOR SLOPE FLOWS 
WITH UNIFORM THICKNESS 

We have obtained a system of equations that is ac- 
ceptable as a model of isotropic pastes. Next, let us 



analyze this system under a particular setup describing 
slope flows. Though the equations are rid-dimensional 
and it is also possible to formulate boundary conditions 
for fully three-dimensional surface deformations, it is not 
wise trying to solve the full system immediately by di- 
rect numerical simulations, as it would require too much 
difficulty and provide too little insight. Rather, the first 
thing to do is to elucidate the basic behavior of the model 
in the simplest situation. 

Here, we study a two-dimensional paste flow (rid = 2) 
on a slope inclined by an angle 9. The setup of the system 
is shown in Fig.E) All of the motion is supposed to occur 
in the {x, z)-plane, and it is in this plane that the paste 
is assumed to be isotropic. The flow is driven by the 
gravitational force F = pG, where G is the gravitational 
acceleration vector. 



G = 



G sin 9 
-Gcos9 



= (Gsm0)e^ - {Gcos9)e^. (33) 



The free surface requires the dynamical boundary con- 
dition and the kinematic boundary condition. The dy- 
namical boundary condition prescribes the continuity of 
the stress, while the kinematic boundary condition postu- 
lates that the surface must move together with the adja- 
cent fluid to satisfy the mass conservation law. Since we 
assume here that the paste layer has a constant thick- 
ness H, the dynamical boundary condition reduces to 
p-'^Vjz = (see Eq. ((X20)l and the text below it in Ap- 
pendix[X|. The kinematic boundary condition is trivially 
satisfied by assuming that the velocity field is also uni- 
form in regard to x and parallel to the cc-axis. 

Under the present assumptions, the fluid motion is ex- 
pressed in terms of a single function which we denote by 
X = X{C,t), as 



r = 



C 



= {£,+X)e.^ + Ce, 



(34) 



The time-derivative of Eq. ((34)) gives the velocity, 

v^dtr = Uea;, (35) 

where U = dtX. By substituting Eq. ([M]) into Eq. (O, 
we obtain the Euclidean metric tensor gij expressed in 
terms of X, 



9ii 5?C 
5C? 5CC 



1 



X' 



X' l + X" 



(36) 



where X' is an abbreviation for dX/dC,. The incompress- 
ibility condition (fTU|) is already satisfied and there is no 
need to require it particularly. 

Let us concretize the equations containing the natural 
metric. The calculation can be performed in at least two 
different ways: one may evaluate the terms in the mo- 
mentum equation (jlip either on the ground of modern 
differential geometry of the Riemannian manifold deter- 
mined by gij , or fully utilizing the Cartesian components 



in the embedding Euclidean space, as is shown in the 
latter half of Appendix [Xj Both methods yield the same 
result. 

As the natural metric tensor for the present case is a 
2x2 symmetric tensor with detg'' fixed to be unity, it 
can be expressed by two parameters. We set 



rtl 



9a 9lc 



(3 






(37) 



with a = a(C,i), /? = /3(C,i), and calculate its inverse 
matrix gt]. Then we substitute it, together with g~^ 
calculated from Eq. ([5S)) . into the equations composing 
the constitutive relation. Eq. P^ then yields the stress 
tensor P. Its Cartesian representation, calculated from 
Eqs. dini) and ^, reads 



P =pl-a 
= pl-S 



S"(l+a2)-l CT 



(38) 



where a — axz/S stands for the nondimensionalized 
shear stress, and is given by 



a = e-"X' - f3. 
The momentum equation (jlip reads 

pdtU ^ sdi:a + pG:, 



(39) 



(40) 



where Gx = G sin 9 is the a;-directional component of the 
gravitational acceleration vector G, given by Eq. ([33)) . 
Note that the depthwise component of the equation of 
motion does not participate in the dynamics, as it deter- 
mines only the hydrostatic pressure. 

From Eq. ([T7| or ^, taking Eq. ([28]) into account 
and using g parametrized as Eq. ([55)1 and g^ as Eq. ([57)1 , 
we obtain 



dta = l- 



TOt 



Tdtf3 =-p + 



2 + e 
2 



2 + e 



-X' 



(41) 
(42) 



where we have utilized the relation K — 2/(2 + e) with e 
defined by Eq. ([^5)). which holds for the two-dimensional 
case (we note that the three-dimensional case is not so 
simple) . By calculating e from Eq. ([^5)) and then rewrit- 
ing the result in terms of a, we find 



£ = e" {e-"X' - pf + 2 (cosh a - 1) 
= e"CT^-|-2(cosha- 1). 



(43) 



Note that Eq. ([43)1 endorses the relation between s and a 
stated several lines before Eq. ([29| . as long as a = o{d) 
(which is usually the case). 

Though the above equations constitute a closed sys- 
tem, X' and (3 are inconvenient variables as they in- 
crease unboundedly as time elapses. To avoid this in- 
convenience, we rewrite the equations in terms of U and 



a. Using the evolution of a instead of Eq. ([42)1 for /3, and 
also rewriting t in terms of v{e), we obtain a system of 
three equations governing three variables, namely a(C, t), 
a(C,i),and(7(C,i): 



d,a ^ He) (l - 1^7 

dtd^e-°'dcU-iy{e)a 

dtU = -dca + Gx. 
P 



(44a) 
(44b) 
(44c) 



Aside from the curious equation (|44ap for a, this system 
of equations has a familiar form that can be recognized 
as a description of a slope fiow (Fig. ^, with Eq. ()44b[) 
relating the nondimensional shear stress a to the shear 
rate d(^U, and Eq. ()44cp describing momentum balance. 
Plasticity is introduced via ^{e) that is the inverse of 
the relaxation time mentioned in Eq. (|T]l . The functional 
form of i'{e) is specified by Eq. (|29p and Fig. [1] on the 
basis of Bingham plasticity. The nondimensional strain 
energy e, defined by Eq. ([^ . is evaluated as a function 
of a and a as in Eq. ([^5)1 . 

Eqs. ([44)) require two boundary conditions. We pose a 
no-slip boundary condition at the wall, 



U\ 



C=o 



0, 



(45) 



while the free-surface condition, for the present case, 
gives 



^\c=H = 0. 



ANALYSIS 



(46) 



A. Qualitative consideration 

Eqs. ([33]) together with two boundary conditions de- 
fine a closed system of evolutional equations. The energy 
is supplied by gravitational work pG^U, stored as elastic 
energy e, and dissipated through the relaxation of g'' that 
represents the viscous part of the Maxwell model. Plas- 
ticity implies that the dissipation process is limited by a 
threshold, in such way that the relaxation time can be- 
come infinitely large according to Eq. ([^5)) . This allows 
some part of the elastic energy to remain frozen inside 
the paste. 

In the present study, a plays an important role. 
Eq. ()44ap clarifies that the threshold mechanism included 
in i^(e), shown in Fig.[l] governs the fundamental behav- 
ior of a. For an e smaller than the threshold value, iy{e) 
vanishes and therefore a practically arbitrary function of 
( is admissible as a steady solution to Eq. ()44ap . as long 
as it allows e to stay within the threshold. This implies a 
strong non-uniqueness of a that can remain in the static 
paste; there are an infinitely large number of possibilities, 
whose realization depends on the time-dependent process 
of evolution (an analogous situation occurs also in dry 
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granular materials subject to static friction [23|)- On the 
other hand, a in the flowing paste is expected to evolve 
toward a steady solution that is uniquely determined if 
the external force, film thickness and paste properties are 
specified. This steady solution will be provided later in 
a closed form. 

We will show that the residence of an a > in the 
paste means the presence of cc-directional tension. Then 
we will derive a steady solution for a flowing paste analyt- 
ically, showing that a is positive there. Time-dependent 
numerical calculations for flowing pastes typically exhibit 
relaxation toward this solution, involving the creation of 
a positive a. The numerical calculations also show that 
some portion of a remains in the paste after its flow is 
stopped, and the residual value of a is still positive. This 
process creates an x-directional tension remaining in the 
paste and therefore gives a possible clarification of the 
Type-I Nakahara effect. 



B. Residual tension 

Let us confirm that a > implies tension. This is 
intuitively evident if we recall that e~" stands for the 
^^-component of the natural metric tensor (j37p . and con- 
ceive of e~" < 1 as contraction of natural length of the 
(supposed) "springs" in the x-direction. More formally, 
this is demonstrated by calculating the normal stress dif- 
ference for the "ground state" that minimizes the elastic 
energy e as a function of g, with g' being fixed. In terms 
of the parametrization given by Eqs. ([55)1 and ([57)l . the 
problem is to minimize e — e{X' , a, (3) for fixed values of 
(a,/3). 

From Eq. ([^5]) we find that the minimizer of e(X', a, (3) 
is X' = e"/3, or equivalently a = Q (vanishing shear 
stress). Then, using Eq. ([55)) to calculate the diagonal 
components of f? in Eq. (|25p . we find the normal stress 
difference 



— 2S'sinhck 



(47) 



for (7 = 0. Clearly, this is positive for a > 0, showing a 
residual tension in the x-direction. 



C. Steady solution for flowing pastes 



an arbitrary a = a{C) such that e < (Ty/S"^ (i.e. i^(e) = 
0). On the other hand, v^e) must be non-zero for flowing 
pastes, which makes the steady solution totally different. 
For steady flows {v{e) ^ and dta = 0), Eq. (|44ap yields 



2(e"-I). 



(49) 



Since e must be positive according to Eq. (j43p . from the 
above equation (|49p it follows that a must be positive as 
well. More concretely, from Eqs. (gS]), (gH]) and (gS]) we 
find 



a 



log (I 
log 



a') 



1 



S 



(H-C) 



(50) 



for the flowing part of the paste in steady state. It is also 
confirmed that e c^ a" for small \a\. 

The neighborhood of the free surface requires a sep- 
arate treatment, because this region remains solidified 
due to the lack of a sufficient shear stress. The boundary 
between the solidified and fluidized regions can be calcu- 
lated by using Eqs. (|i51) and ([SD]) . which give e = e{a{()) 
in the fluidized region, to find the location Cy such that 
e(a(CY)) = (Jy/S'^. In the region (y < ( < H where the 
paste is solidified, the velocity U is uniform. The velocity 
U in the fiuidized region can be obtained by integrating 
Eq. (j44bp under the boundary condition ((45|) . 

As is evident from Eq. (|48p , the maximum of the shear 
stress (Txz occurs at the wall. The wall shear stress and 
its nondimensionalized value are (Txz\r^Q — pGxH and 
^lr=o ~ pGxH/S. For the paste to flow steadily, this 
wall shear stress pG^H must be greater than the yield 
stress ay. The maximum a also occurs at the wall: 



1 



max a 

C 



log 



pGxH 

S 



(51) 



according to Eq. ([50)) . 

The above discussion suggests two nondimensional pa- 
rameters that can be expressed as a ratio pGxH : ay : S 
Let us complete the dimensional analysis of Eqs. ([44] 
before proceeding to the numerical calculation of time- 
dependent solutions. 



Eq. ([47)1 shows that a paste layer left in the unloaded 
state ((7 = 0) bears an x-directional tension if a > 0. 
The next task is to show that the flow makes a > if it 
approaches a steady solution of Eqs. ([1^ . 

For steady flows, the nondimensional shear stress a — 
f^xz/S is determined by the momentum balance (|44cp and 
the free surface boundary condition (j46p . The result is 



pGx 



(H-C). 



(48) 



Note that Eq. (|48p holds for static states as well. For that 
case, the steady solution consists of Eq. (HS)) . C/ = 0, and 



D. Dimensional analysis 

Eqs. ()44p contain five physical parameters, namely S, 
rjp, p, Gx and cry (the last one comes through v). The 
first three determine the viscoclastic time scale tq = Vvl ^ 
and the length scale Iq — rjp/^/pS. The boundary condi- 
tions introduce the layer thickness H as another length 
scale. 

The system is characterized by three nondimensional 
parameters, for example, H/£q, ay/S, and pGxH/S (or 
a suitable combination of them). Note that pGxH gives 



an estimation of the wall shear stress, which represents 
the magnitude of the external forcing. 

Evaluation of Reynolds number will be useful for con- 
sidering the Newtonian limit. On the basis of rjp, H, and 
U ~ pGxH/{rip/H) — pGxH^/rjp, it is estimated as 



R 



H ~ ^ 



pGxH 

S ' 



this is indeed calculated from two of the three parameters 
stated above. In the present setup, H is taken as the 
representative length scale, but we point out a general 
possibility that the system may be characterized by other 
Reynolds numbers, such as Rl = UL/rjp based on the 
horizontal length scale L. In future studies this point 
may have to be taken into account. 



E. Numerical calculation of unsteady solution 

The author calculated the numerical solutions of 
Eqs. p4)) (slightly modified, as we will see below) un- 
der the initial condition (a, a, U)\^^q — (0, 0, 0) and the 
boundary conditions (|^5|l and (|46p . with iy{e{a,a)) de- 
fined by Eqs. (P5|) and (P5|) . for hundreds of different 
nondimensional parameters. With the hyperbolic char- 
acter of Eqs. (HH) taken into account, the calculation 
adopted the two-step Lax-Wendroff scheme [2J| . 

Since we are interested not only in the creation process 
of a, but also the storage of a after the flow is stopped, 
it is necessary to simulate the process to stop the flow. 
To this aim, we "switch off" gravity at some time t ~ 
T* (» To), replacing Eq. ((iic)) by 



dtU 



S ^ . Fx 
-dra + — 

P P 



(0 < i < T*) 

(t > n) 



(52) 



with T, = IOOto or T* = 200ro. 

Fig. [3] depicts a typical evolution of (a, ct, f7). The pa- 
rameters are H — 5£o, pGxH : ay ■ S — 0.6 : 0.3 : 1, 
and r* — IOOtq. In the first stage of the evolution, the 
system rapidly approaches steady state, except for the 
region adjacent to the boundary between the fluidized 
and solidified regions {( — (y = 2.58 io) where the re- 
laxation time is significantly longer. After the gravity is 
"switched off" at t = T* , both a and U oscillates around 
zero. This oscillation should be damped if we consider 
the solvent viscosity, which is neglected in the present 
model. What must be noted is that a remains finite, 
though it decreases, after the driving force is switched 
off at i = T*. The sign of the residual a is positive. 

According to the analysis of 801 cases with T = IOOtq 
and 689 cases with T — 200to, the behavior of a for differ- 
ent values oi pGxH is summarized as follows. For pGxH 
smaller than cty/2, throughout the evolution a remains 
zero. If pGxH exceeds crY/2 but still remains below cry, 
the evolution during the forcing (0 < t < T*) is basically 
unsteady, where a is produced little by little from the 



interference of the stress waves. For ay < pGxH < S 
(we always assume ay < 5), steady yield flow occurs, 
creating a according to Eq. ([50)1 . In both regimes stated 
above, a residual a is observed after t — T^,. The steady 
solution, Eq. ((50)) . ceases to exist for pGxH > S, which 
leads to the unlimited acceleration of the flow. This last 
case is out of the scope of the present model, because U 
should be limited if, again, the solvent viscosity is taken 
into account. 

Steady solutions obtained by the time-dependent cal- 
culation during the forcing, approximately for T^/2 < 
t < (3/4)r*, are checked against the analytical solution 
in Fig. U The curve shows max^ a given by Eq. (j51[) as 
a function of pGxH/S. The symbols, consisting of 122 
circles (T* = lOOro) and 110 triangles (T* = 200ro), in- 
dicate the numerical values of max^ a calculated within 
the range ay < pGxH < 0.72 5", 0.05 5" < ay < 0.30 5, 
and 0.20 4 < H < 8.0 4. The size (and the color) of 
each symbol indicates the magnitude of ay/S. Fig. [3] 
demonstrates that the value of steady a is independent 
of ay, once pGxH exceeds it. For T^, = IOOtq, there were 
several cases for which a did not attain its steady value 
(with the criterion Aa/|a| = 0.02); these cases are elim- 
inated from Fig. [4] for clarity. Such exceptional cases did 
not occur for T* = 200to. 

As an explanation of the Nakahara effect, it is essential 
to show that some of a remains in the paste even after 
the flow is stopped, instead of decaying away. Fig. [5] 
shows the numerical values of maxij a remaining steady 
(not to decay any more) after the ffow is stopped. Here 
not a itself, but aS'^/ay is plotted against pGxH/ay for 
pGxH > ay/2 (the ranges of ay/S is the same as in 
Fig. H and that of H/Iq is 0.05 4 < H < 8.0 Iq). The 
values of residual a for ay < pGxH < 1.3 ay is fltted by 



max a 

C 



S^ 



0.75^^^-0.2 
ay 



(53) 



In contrast to the steady value of a in the flow subject 
to the driving force, the residual value of a in Eq. ((53)) is 
strongly dependent on ay. In particular, if pGxH/ay is 
kept constant, Eq. (|55|) states that the residual value of a 
is scaled by {ay/S)^. This result seems understandable 
if we assume that, during the decay of a and a = axz/S, 
the first equal sign in Eq. ([50]) remains valid, until axz 
reaches the threshold value ay. This gives a rough esti- 
mation of the residual a ~ 0.5(o'y/5)^. Unfortunately, 
theoretical clarification of Eq. ([53]) in regard to its de- 
pendence on pGxH/ay is not currently available. 



VI. DISCUSSION AND CONCLUDING 
REMARKS 

A. Relationship with crack pattern experiments 

In this paper we have found the creation and fixation of 
the x-directional tension using a model equation for flows 
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FIG. 3: Typical evolution of {a, a, U). 





pGxH/oY 



FIG. 4; Steady values of a during the forcing. Spatial max- 
imum of a, whose steadiness is checked for a certain time 
range, is plotted against pGxH/S. The line shows the ana- 
lytical solution (|5ip . The circles and the triangles represent 
numerical values for T* = lOOro and T, = 200to, respectively. 
The size and the color of the symbols indicate uy/S, from the 
small blue symbols for cty = 0.05 5* to the large green symbols 
for CTY = 0.30 S. 



FIG. 5: Residual value of a that remains in the paste after 
the forcing is removed. The spatial maximum of a, rescaled 
by (cty/S)^, is plotted against pGxH/aY- The same symbols 
(circles and triangles) are used as in Fig. 3] The nonverti- 
cal broken line, with slope 0.75, represents the fitting rela- 
tion ((55)) . 



of isotropic pastes. This provides a possible scenario for 
the Nakahara effect (Type I). 

During the drying process, the paste slowly shrinks. 
Mathematically, this process is described as an isotropic 
contraction (shrinking) of the natural metric g^. If the 
paste had not undergone a flowing process, this contrac- 
tion would produce a basically isotropic tension in the 
(x, y)-plane (parallel to the surface and the bottom) and 
therefore would lead to isotropic crack patterns. Actu- 
ally, this is not the case: we have found that a positive a 



is created during the flowing process, which implies that 
the natural metric is already contracted in the stream- 
wise direction. Strictly speaking, the present analysis is 
limited to the two-dinrensional system in the (x, z)-plane 
and therefore it cannot tell whether any y-directional 
contraction occurs, but it is unlikely that it will occur 
to the same extent as the a;-directional contraction. In 
fact, though a full analysis of three-dimensional system 
is too complicated to develop here, a simple perturba- 
tion analysis supports the above conjecture. The bonds 
perpendicular to the flow are therefore the first ones to 
break, causing cracks perpendicular to the flow and thus 
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clarifying the Nakahara effect. 

The present numerical analysis predicts that the mag- 
nitude of the residual a is scaled by {ay/S)^, as is seen 



in Fig. [5] and Eq. ((53)1 . This result is consistent with the 
observation of Nakahara and Matsuo in regard to the 
strength of the memory effect summarized as Fig. 2 in 
Ref. Q. The figure presents a classification of the ob- 
served patterns as a function of the solid volume fraction 
(density of the paste) and the strength of the external 
forcing. Its Region B, which lies just above the yield 
stress line and exhibits the memory effect, is subdivided 
according to the strength of the anisotropy in the pattern; 
strong anisotropy is observed for denser pastes (lamel- 
lar crack patterns, denoted by solid squares, occupy the 
subregion with volume fraction greater than 40%), while 
less dense pastes exhibit weaker anisotropy, resulting in 
large-scale lamellar cracks {':^ H) combined with cellular 
structure with smaller length scales (~ H). If we ad- 
mit that cry/S is greater for denser pastes, the difference 
in the strength of anisotropy can be explained from our 
theory predicting a oc {uy I SY. 



B. Comparison with dry granular flows and other 
systems exhibiting memory effects 

Memory effects are quite common in many glassy sys- 
tems, ranging from granular matters to spin glasses. 
In the case of dry granular matters [23, [23], history- 
dependent behavior essentially originates from the ex- 
istence of interparticulate static friction. According to 
Coulomb's friction law, the interparticulate forces admit 
static indeterminacy, giving rise to the history-dependent 
stress state. Fluidization and solidification of granu- 
lar matter also involves the creation and destruction of 
grain-scale structures, such as arching and force chains. 
Though the present study on paste flows is based on 
macroscopic description and therefore discussion on the 
grain-scale structure is outside its scope, comparative 
consideration on static indeterminacy is quite helpful 
in understanding some common mechanisms underlying 
paste flows and dry granular flows. 

As seen at the top of Sec. |Vl the threshold mechanism 
in v{e) results in the static indeterminacy of a. It is 
this static indeterminacy that enables the retention of 
the memory of the shear flow. (In a different setup [l3| , 
residual stress is introduced via static indeterminacy of 
/3.) Thus the present paste model shares an important 
feature with dry granular systems. 

To elucidate the analogy and distinction between the 
Bingham plasticity and Coulomb friction, let us consider 
an instructive problem taken from Chapter 3 of Duran's 
book [23]. Suppose a brick on an inclined wall, subject to 
static Coulomb friction (coefficient /Xg) and a spring, as is 
illustrated in Fig.[6Ka). Duran's problem is to determine 
the deformation x (or equivalently the repulsion kx) of 
the spring as a function of the inclination angle 0, when 
varies slowly in time. 



Suppose that the wall starts from the horizontal po- 
sition {9 = 0) and that we know the initial value of x, 
which we denote by xq . For a while x is stuck to xq , until 
the "yield" criterion 

\mG sin — kx\ — ^s'niG cos 9 

is attained and the brick starts to slip. We assume the 
viscous resistance —ex and neglect the dynamic Coulomb 
friction for simplicity [50| , so that the brick moves accord- 
ing to 

d^a; dx ^ ■ /, n • i • • \ 

m — 7T = —c— kx + mG sm brick m motion 

dt^ dt 

and eventually stops. This process is repeated while 
increases, as is shown in Fig. [Tfa) with a solid line (each 
slip is assumed to stop when kx — mG sin according to 
Duran [23]). If 9 starts from 7r/2 and decreases slowly in 
time, a similar but different stick-slip motion occurs, as 
depicted by the broken line. Thus, the system exhibits 
mechanical hysteresis due to static friction. 

Now let us compare this mechanical hysteresis with the 
behavior of the system in Fig. (Sfb) , where the Coulomb 
friction is replaced by a discrete-element analogue of 
Bingham-like elastoplasticity. Its behavior is defined by 
combining Eqs. ([H]) and (jlSp with 



KT) 



(]£^(r)] < threshold) 

T^^ (]£;(r)] » threshold) 



(54) 



which is a discrete-element version of Eq. ([T]), with 
E{T) = T^/{2k) standing for the elastic energy stored 
in this element. The governing equation of this system is 
summarized as 



dr^ 
^^,iT){x~x^) 



m- 



d^ 
di^ 



-T — kx + mG sin ( 



(55a) 
(55b) 



supplemented with T = k{x — a;^). As for z^(T), a 
Lipschitz-continuous form analogous to Eq. ([29]) is as- 
sumed. By numerical integration of Eqs. (j55p with 
increased slowly from zero to 7r/2 and then decreased 
back, we obtain the result shown in Fig.[7]Jb). The thick 
solid line indicates that a shift of a;^ has occurred during 
the process of increasing 9, and this shift was not re- 
covered at all when 9 was decreased (thick broken line). 
In addition, when 9 has returned to zero, there remains 
a difference in x and x\ indicating residual pressure in 
this case. Thus, again, a hysteresis due to static inde- 
terminacy is observed. There is an important difference, 
however, that the curves in Fig. Wijo) are much less sin- 
gular than those in Fig. (Tjja). In other words, at least 
for the values of the parameters and the functional form 
of i^(T) used in this calculation, no stick-slip behavior is 
observed. This is probably related to the property of the 
Bingham model, which predicts continuous shear stress 
across the yield front in quite general cases [2a ]. 
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It is an interesting attempt to reformulate the stick- 
slip motion subject to static Coulomb friction in terms 
of v, to obtain a (formally) unified equation: 



^ = K.,F,iV)(^_. 



(56) 



with 



da; 



V — ——, F — —kx + mG sin 9, N = mG cos I 
at 



A naive choice for v is 



iy{v,F,N) 



{v = and |F| < ^i^N) 
c/m (otherwise). 



Though this function is too singular to constitute a math- 
ematically sound evolutional equation, adoption of con- 
tinuous interpolation similar to Eq. (|29p enables the nu- 
merical integration of Eq. (|56p , resulting in stick-slip mo- 
tion shown in Fig.[7Kc). Note that, in Eq. ((56|) . the relax- 
ation is attributed to the momentum, but this seems to 
be somewhat unnatural if we consider that friction is a 
property of the interface while momentum concerns the 
whole mass of the body. Rather, in analogy to Eqs. ([55]) 
where relaxation is attributed to x\ it seems more ap- 
propriate to introduce a variable describing the state of 
the interface (possibly similar to the one introduced by 
Carlson and Batista [23) and prescribe its relaxation. 
This is beyond the scope of the present work, however. 

In soil mechanics, a continuum version of Coulomb fric- 
tion is known as Mohr-Coulomb plasticity [28] . Its appli- 
cation to the statics of granular materials is usually sup- 
plemented with the limit-state assumption, which states 
that the ratio of the shear stress to the normal stress is 
just below the threshold value everywhere. This assump- 
tion makes it possible to evaluate the stress field with- 
out introducing granular elasticity that is not understood 
very well. However, this theory encounters a number of 
difficulties, as is discussed by Kamrim and Bazant [23 ]. 
According to this theory, the static stress field is sub- 
ject to a nonlinear hyperbolic system of equations (not 
in space-time but in the (x, y)-plane), which predicts a 
highly discontinuous stress field. The solution for the 
velocity field can be even more singular, which seems ab- 
normal both physically and mathematically. Sometimes 
it also fails to satisfy the boundary conditions. Kam- 
rim and Bazant [23] have shown that these difficulties 
can be avoided, within the framework of Mohr-Coulomb 
plasticity with the limit-state assumption, by introducing 
diffusive motions via mesoscale objects called "spots." In 
spite of this successful result, the theory is not free from 
the limitation due to the limit-state assumption, as the 
authors themselves admit that clearly it breaks down in 
some cases. 

Kamrim and Bazant [29!] state repeatedly that the in- 
troduction of elasticity will solve the difficulties of the 
Mohr-Coulomb plastic model. To some extent, this re- 
mark applies to Bingham plasticity as well. For example, 



the original Bingham model exhibits a singular behavior 
due to the lack of elasticity, in the sense that the propa- 
gation speed of yield front is infinitely large [30] • Treat- 
ment of residual stress would be also very difficult, if not 
impossible, without considering finite elasticity. This is 
why we sought to develop an elastoplastic paste model 
from the beginning. 

The present theory is conceptually akin to the mod- 
els of the memory effect in polymeric materials [1^ [31| . 
Miyamoto et al. [l3] studied the memory effect in the 
glass transition of vulcanized rubber. They explained 
their experimental results with a Maxwell-like model. 



O-(i) = a-ruhhcr{'T{t),j{t)) 

ft 

' '-'glass 



bit)-lit')]^dt', (57) 



where a is stress, T is temperature, 7 is strain, Q( 
normalized relaxation function, and t, defined by 



IS 



i=i{t,t')= [ 



du 



T{T{u),-i{u)) 



stands for the intrinsic time lapse. The effect of tem- 
perature control (quenching and reheating) is expressed 
via r, which changes the pace of the intrinsic time t' 
and thereby affects the memory function in Eq. (j57|) . 
Note that 7(t') in the integral can be read as the nat- 
ural length of a spring born at the time t' . In this 
sense, Eq. (|14|) can be regarded as a simplified version 
of Eq. ([57]) , though there is an important difference that 
the memory in Eq. ([TJ]) is ascribed to a single variable 
x^ , while Eq. (|57p can memorize more about the history 
of 7(i'). The "memory capacity" of Eq. ([57]) depends on 
the property of the relaxation function G{ - )■ Using a 
sum of two exponential functions, which implies double 
relaxation, Miyamoto et al. [lO] has successfully repro- 
duced the memory effect, including the effect of aging. 
In the present model, contrastively, the relaxation time 
T is assumed to be a single scalar function. Instead, the 
spatial distribution of a and the effects of nonlinear elas- 
ticity are taken into account, thus enabling the creation 
and storage of the streamwise tension. 

The idea of ascribing the memory to the plastic shift 
in the neutral point of elasticity, corresponding to a in 
Eqs. ((2S1) and (|iia|) . x^ in Eqs. ^ and ((55a|) . and -f{t') 
in Eq. ([57]) . is also shared by Ohzono et al. [31[. They 
studied microwrinkle patterns produced on a platinum- 
coated elastomer surface, governed by the competition 
between the restoring force of the platinum tending to be 
less curved and that of the elastomer that aims to shrink 
back. At room temperature, application of a uniaxial 
compression force breaks the force balance and changes 
the wrinkle pattern, but the original pattern is retrieved 
after the external force is removed. Contrastively, a pro- 
tocol involving higher temperature (annealing-cooling- 
unloading protocol) changes the wrinkle pattern, intro- 
ducing strong anisotropy. The new pattern is less sta- 
ble to external forcing at room temperature, suggesting 
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Bingham-like 
element 




Cu.iiu,, 

frictioi 



FIG. 6: Duran's brick on a slope to illustrate static inde- 
terminacy, (a) Original setup due to Duran 23]. Besides 
the gravity and the repulsion of the spring, the brick is sub- 
ject to static Coulomb friction, (b) A modified setup, where 
Coulomb friction is replaced by a Bingham-like elastoplastic 
element. 



the presence of multiple metastable states. These exper- 
imental results are compared with a model prescribing 
the minimization of the elastic energy as a functional of 
the surface elevation z — z{x,y), 

^ [^\ ^bending ' ^in-plane i ^substrate 5 

C^substratc = / {a {z - Zm)^ + b {z - Zy^)*} dxdy , (58) 

where C/bcnding and C/in-piano are the potentials of bending 
and in-plane deformation of the platinum layer, C/substrato 
is the potential of the substrate (with the constants a and 
b specified explicitly in terms of the material constants of 
the elastomer), and Zm = Zm(a;, y) represents the neutral 
point of C/gubstratc- Correspondence with Eq. (|14p is obvi- 
ous. The memory is carried by spatial distribution of Zm, 
which is fixed at the room temperature but is subject to 
plastic flow in the annealing-cooling-unloading protocol. 
Generally, in elastic systems with more than several 
degrees of freedom (and particularly in continua), a shift 
in the neutral point introduces mechanical frustration. 
In the case of Ohzono et al. 31], it modifies the existing 
frustration, introducing multiple stability. In addition, 
spatial heterogeneity of a and /3 in Eq. (|57|) is equiva- 
lent to the continuous distribution of edge dislocations 
and screw dislocations, respectively [ij, [S^, [33|. Thus, 
frustration is observed universally in systems admitting 
plasticity (in any sense of the word), ranging from gran- 
ular matters to metal crystals and spin glasses. From 
this viewpoint, we understand Eq. (j44ap as describing 
the dynamic creation and static retention of mechanical 
frustration, presenting a macroscopic analogue of dislo- 
cation dynamics. 



C. Future directions 

The present study is entirely based on macroscopic 
phenomenology. It predicts the presence of a macro- 
scopic mechanism that leads to the Type-I Nakahara ef- 
fect, but it does not assert the absence of other mecha- 
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FIG. 7: Mechanical hysteresis in Duran's brick. The deforma- 
tion X is plotted against the inclination 9, with a solid line for 
the forward process (increasing 6) and with a broken line for 
the backward process, (a) Duran's solution [23| in the case of 
static Coulomb friction, (b) Bingham-like case with Eqs. (|55|l . 
In addition to x, the natural length of the Bingham-like ele- 
ment, a; , is delineated in thick (red) lines, (c) A solution of 
Eq. (|56|) exhibiting Coulomb-like stick-slip motion. 



nisms, such as the creation of bond fabric or microscopic 
texture. A possible scenario is that the microscopic bond 
structure is well represented by the macroscopic (hydro- 
dynamic) variables, such as a and a, so that the most 
important feature of the mechanism is already captured 
by the hydrodynamic equations. In other words, we ex- 
pect something analogous to ferromagnetism, where the 
macroscopic magnetization represents the order param- 
eter. We cannot deny the possibility that some pastes 
have "anti-ferromagnetic" bond structure, which makes 
the macroscopic description more difficult. Even in the 
"ferromagnetic" case, consideration of microscopic de- 
tails may introduce some modification. For example, we 
have regarded the yield stress as a given constant, but 
this may possibly need to be modified, in the way similar 
to work hardening and the Bauschinger effect in met- 
als [34!| . The constitutive relations assumed in this paper 
require justification more certain than a physicist's in- 
tuition, either by microscopic analysis or thermodynam- 
ical inspection. It is worthwhile to consider an exten- 
sion of microscopic theories for glassy liquids, such as 
the mode-coupling theory [33, [3a| and the pair distribu- 
tion function theory [33], in the direction corresponding 
to that of direct-interaction approximation in fluid tur- 
bulence based on Lagrangian description [3a, |33, EO, Ell , 
in search of microscopic expression for g^J . Such a mi- 
croscopic approach would also allow us to construct a 
model for pastes whose properties are not isotropic. It is 
expected, for example, that a model including competi- 
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tive interaction between the natural metric and director 
field may clarify the Type-II Nakahara effect. 

Within the framework of the present model, an ex- 
planation of Eq. ((53)) that gives the magnitude of the 
residual a is an open question. It is also necessary to ex- 
tend the present work in several respects. On one hand, 
the limitation to uniform flows must be removed. The 
stopping process is simulated in this paper by switching 
gravity off, but in real experiments the paste flow stops 
when the paste supply is cut. Simulation of this pro- 
cess requires the introduction of a variable layer thick- 
ness h — h{x,t), where the flow and the stress fields 
depend on x as well. This extension will clarify the rele- 
vance of different mechanisms, such as the one proposed 
by Otsuki 42| where the a;-dependence of the plastic de- 
formation is essential. The present model can be readily 
extended in this direction, though its numerical analysis 
will be much more difficult. Derivation of reduced equa- 
tions, such as depth-averaging (corresponding to Shkadov 
model [4a . |44| . |45| in the film fiows and Saint- Venant 
model [46l . |47[ in civil engineering) , will be worth consid- 
ering. 

On the side of experiments, it is desirable to realize a 
uniform slope flow by eliminating the boundary effect in 
the y direction. It is also necessary to measure the paste 
properties, such as S and ay, so that a qualitative com- 
parison between the theory and the experiment becomes 
possible. Finally, since the mechanism proposed in this 
paper is closely related to the nonlinear viscoelasticity, it 
will be highly supportive to detect any indication of non- 
linear viscoelasticity in the paste, such as Weissenberg 
effect. 
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APPENDIX A: NOTES ON THE FORMULATION 

OF LAGRANGIAN CONTINUUM MECHANICS 

IN TERMS OF DIFFERENTIAL GEOMETRY 

Here we summarize the minimal mathematical knowl- 
edge required to understand, for example, how to cal- 
culate each side of Eq. pT|) . Instead of going along the 
rather expensive highway of Riemannian differential ge- 



ometry, we take a shortcut, making full use of the n^- 
dimensional Euclidean space where the whole system is 
embedded. 



Contravariant vector components 

For each instant (with t fixed arbitrarily) , the mapping 
from ^ to r provides an instantaneous curvilinear coordi- 
nate system. This is sometimes rcfcrcd to as a convected 
coordinate system [48| . It is this coordinated system, and 
not the space, that is curved. 

Provided that the mapping (O, for fixed t, is suffi- 
ciently smooth and locally invertible, we find that 



dr dr dr 

9^ ' drj' d( 



(for rid — 3) 



forms a set of local bases in the nd-dimensional Euclidean 
space (r-space). Then, an arbitrary vector field, say f, 
can be expressed as 



f = 



dr dr dr 
d(, drj d( 



/" 






or, in abbreviation with Einstein's contraction rule, 

f = Pd.r. (Al) 

The coefficients (/*) in Eq. (jAip are referred to as con- 
travariant components of the vector field f . According to 
the convention of differential geometry, the contravariant 
components are superscripted. 

If the labeling variable is changed from ^ to ^ (in terms 
of a continuous, one-to-one mapping independent of i), 
the bases are changed to 



dr 



d^J dr 
'd^W 



Meanwhile the change from (/') to (/*) occurs in such a 
way that it cancels the change in the bases (therefore the 
name "contravariant" ) , so that the vector f itself remains 
unaffected: 

f = p^= p^ 

■I g^i J d^t- 

An equation describing the relations between physical 
quantities should be independent of the choice of label- 
ing variables. This is assured if and only if every term 
on both sides of the equation has the same behavior in 
regard to the relabeling. For example. 



= 26* 



is acceptable, while 
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is not (we cannot add a scalar 2 to a contravariant vector 
component V). 

A second-order tensor, say P, can be expressed as 

P = P'^ {d.,r) (g) (djr) (A2) 

where (g) is the tensor product, such that 

a • (b (Ki c) = (c (g) b) • a = (a • b) c. 

The contravariant components (P'^) are subject to the 
same kind of change as the product of two contravari- 
ant vector components, so that P remains unaffected by 
the relabehng. Note that Kronecker's delta with super- 
scripts, S'^ , does not behave properly in regard to re- 
labeling and therefore is not acceptable as a physically 
meaningful tensor. 

Dual basis and covariant vector components 

As has been stated, we assume that the mapping from 
^ to r is smooth and invertible. Therefore, it makes sense 
to define 



or 



(A3) 



a nabla without subscript, V, is a mere abbreviation for 
d/dr, i.e. the gradient operator in the r-space. Evidently 
{V^'} is the dual basis of {div}: 



(5,r) • Ve = S,^ 
due to the chain rule. Also 

(Vf ) (g) 9,r = 1 



(A4) 



(A5) 



where 1 denotes the unit tensor in the r-space. Note that 
Eq. (jASP holds thanks to the fact that the embedding r- 
space has the same dimension as the ^-space (otherwise 
(V^') (g) diT would be a projection operator whose rank is 
lower than the dimension of the r-space) . The dual basis 
allows us to find the contravariant components of a given 
vector field, say f , by 



r-f-vf; 



(A6) 



substitution of this expression into the right-hand side of 
Eq. ((XT|) recovers f due to Eq. ([Xsjl . 

As opposed to the contravariant components (/') of a 
vector field f , we define its covariant components (/i) by 

f = /.Vr. (A7) 

It is easily confirmed that 

/, = {d,r) ■ f = g,,f 

where gij stands for the Euclidean metric tensor defined 
in Eq. ©. 



Using Eqs. ([8]) and (jASp . we identify (gij) with covari- 
ant components of the Euclidean unit tensor, 

g^J (vr) ® (ve) = a. 

The contravariant components of 1 comprise the inverse 
matrix of (gij), denoted by (g'^), which leads to Eq. (|24p . 
This imphes that (pg^^) in Eq. ^ and {Kg'^) in Eq. HH) 
stand for isotropic tensors. 



Covariant derivative 

The momentum equation (jlip , represented in terms of 
contravariant components, contains Vj which generally 
differs from dj = d/d£,^ . This "nabla with a subscript" is 
referred to as the covariant derivative. When the space 
is curved, it is not a trivial problem to define the covari- 
ant derivative in an appropriate way. Fortunately, since 
the space itself is now flat, we can now define Vj as a 
component of a simple "gradient" using V = d/dr. For 
a scalar field, say (p, its gradient is 



V(p = 



dip 9^' dip 



dr 



dr d^^ 



i^Od^ip; 



the covariant derivative of ip is given by the covariant 
components (i.e. the coefficients for V^*) of V(p, 



V^ip = diip. 



(A8) 



The covariant derivative of a vector field is slightly more 
complicated. For f given in terms of its contravariant 
components (/*), the gradient is 

gradf = V (g f = ((V^^) a^) (g ifd^r) 
= (Ve-'')®9,(r9.r); 

we define the covariant derivative Vj/* by 

d,ird,r) = iV,r)d,r (A9) 

so that 

gradf =(V,r)((Ve)®9,r). 

A handy way to evaluate ^jf^, in the present case, is 
to calculate the Cartesian components of f = Pdir and 
then to differentiate them with ^^ , which yields the left- 
hand side of Eq. (|A9p . For those who disdain to depend 
on the embedding r-space, there is a more orthodox way 
based on a formula 



y^^ = dJ' 



^L/^ 



with F;^^ referred to as Levi-Civita connection (also 
known as Christojfel symbol when it is calculated from 
didjr). Both ways lead to the same result. 
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The momentum equation (fTTj) contains a term arising 
from the divergence of stress tensor, 



div P — Hm ■ 



AV 



P-ndS-^^ V-*P 



d{AV) 



where *( • ) denotes transposition (practically it could 
be omitted, as P is symmetric). Substitution of P = 
P'^d^r) ® (djr) and V = {V^'')dk yields 

d 



divP = — •t(P'^(ft;r)®(9,r)) 



(AlO) 



where the last equal sign follows from the definition of 
the transposition. At this stage, we need the covariant 
derivative for (P'-'). Taking into account a general postu- 
lation that any formula for a second-order tensor should 
apply to the tensor product of two vectors as well, we 
find the appropriate definition to be 

dk (P'^d^r) ® (a,r)) = {\/kP'') (5,r) ® (d^r) (All) 

so that 

div P = (VkP'n mr) (g> (djr)) ■ (V^*'-) 
= (VfeP^^)(9,r)5/ 
= (V,P'^)a,r. 

Again, VfcP'^ can be evaluated either in terms of the 
Cartesian components of P or with a formula 



VfcP^-' = SfcP'- 



niP'' 



niP" 



unless specified otherwise (in Eq. ([4]), for example). The 
velocity v is then given by Eq. ([6|), and the (material) 
acceleration is 



d^r = dtv = dt (y'^^^r) 



(A12) 



as is seen on the left-hand side of the momentum equation 
just above Eq. ((TT|) . Taking the time-dependence of diY 
into account, we evaluate the acceleration as 

= {dtv'')d^v + v''d^dtr 

and rewrite the last term, which contains 9jV, with the 
covariant derivative. Thus we find 



dt-v = {dtv' + vi\/jv') d,v. 



(A13) 



The contravariant component of Eq. (jA13|l . multiplied 
by p, gives the left-hand side of Eq. pTjl . 



Derivation of Eq. (fJO]) 

Next, we study a concrete example to see how the mo- 
mentum equation (jlip is evaluated. With the mapping 
^ 1^ r specified as Eq. (f34|) , the momentum equation pTJ) 
is to be reduced to Eq. (|^0)l . 

Eq. (|34p readily yields the velocity in Eq. ([35|) and the 
local basis 



Velocity and acceleration 

Up to the present point in this appendix, we have 
treated the spatial aspect of the mapping from (^, t) to 
r with t fixed. Now we will detail the temporal aspect 
of this mapping. Let us recall that dt stands for the La- 
grange derivative, 



dt{-) 



d_ 
'dt 



d(V 



dcv 



X' 

1 



where U and X' are understood as 



U^dtX{Q,t), x'^d^xict). 



(A14) 



Substituting Eq. (|Al4j) into Eq. ^ yields (g^) in 
Eq. 



The natural metric tensor is parametrized as Eq. p7p ; this expression becomes identical to that for g if a = and 
/3 = X' . The components of the inverse natural metric tensor are then 



4' 4' 



{l + (3^)e°' -13' 
-13 e-" 



(A15) 
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By using Eqs. (|A14p and (jAlSp . the term gl^div) (g) (djr) in Eq. ([25]) is calculated to be 

gl'id^r) (g (djr) = (1 + /32)e"(9^r) ® (^^r) - /3 ((a^r) ® (dcr) + (dcr) ® (^^r)) + e-"(9cr) (9cr) 



'2\„a 



il + P')e 



1 




p 



2X' i 
1 



X'^ X' 
X' 1 






(A16) 



r 



Taking notice of the (x, z)-coniponent of this expression, 
which corresponds to Uxz/S, we introduce a given by 
Eq. (|55)l . Then Eq. (P5|) yields a concrete expression for 
CT shown in Eq. ([38]). 
In the present setup, the nabla operator is given by 



^ = IT = ^^? 
or or 



ac 

5r 



ar 



(A17) 



where 






1 



ac 



Then the divergence in the momentum equation (jlip is 
evaluated in terms of the Cartesian components in the 
r-space: 



- div P = - 



1 



d. 



'iP^ 



a, 



c 
a 



'C-P 
1 CT 



1 
-X' 1 






s 






where it is taken into account that a, j3 and X' are inde- 
pendent of ^. As for the left-hand side of the momentum 
equation, it is easily shown that 



Stv 



dtU 




Calculating the inner product of the momentum equa- 
tion with a^r, we obtain 

pdtU = -a^p + Sd^a + pG sin 6; (A18) 

similarly, the inner product with Oqy yields 

pX'dtU = -dc_v + S (X'dQd- - e^^aca) 

-t-/9G'(X'sin6'-cos6'). (A19) 



where Patm denotes the (constant) atmospheric pressure, 
which can be set equal to zero without loss of generality, 
and rij stands for the covariant component of the surface 
normal vector, which is given by n = V(z — H) so that 
Uj = VjZ = dz/dS,^ for the present case. For p, the 
boundary condition (|A20p reads 



{p-(^zz)\ 



C=H 



Patm (= 0) 



(A21) 



with ffzz — S{e~°' — 1) according to Eq. (|55|l. Evidently, 
Eq. (|A21|) is also independent of ^. Then p turns out 
to be totally independent of ^, which implies that d^p in 
Eq. (JAISP vanishes, leading to Eq. (|i(I)) . 



Derivation of Eqs. (gT]) and (|42|) 

The relaxation of g'' is described by Eq. (fT7|) or, equiva- 
lently, Eq. P^ . We substitute g parametrized as Eq, 
and g^ as Eq. ^T^j into Eq. (HH), together with 



Tdt 






(l+/32)e° 




rata 



2/3e" -1 
-1 



Tdtp. 



Equating each component of the matrix yields three equa- 
tions for two variables a and /3; the equations are consis- 
tent (solvable) only when K is set appropriately, which 
is calculated, according to Eq. ((28|) . as 



K 



2 cosh ( 



(A22) 



with e given by Eq. (j43p in the two-dimensional 
case. From the ('('-component and the C^-component of 
Eq. IH]) we obtain Eq. (gT]) and Eq. ^, respectively. 



From Eq. (|A19p . we find that d^p is independent of ^. 

Here we use a concrete formulation of the free-surface 
boundary condition for P (neglecting surface tension and 
surface contamination) , 



P'^n 



'1\C=H 



Patmff^-'nj 



(A20) 
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APPENDIX B: VARIATION OF THE ELASTIC and 

ENERGY E 



Eq. (|T6|) is obtained from elastic energy E in Eq. ((23|) 
by calculating its variation in regard to r = r(^) under 
the constraint det g - 1 . In this calculation we use q^ y^i^ ^ i y^^ ^ki q^^^^ ^ y^^ ^y^k ^ . q^ q^^_ 



2 



(5(detg) = (detg)5*^%j, 
The result is as follows: 



and 



= -S j ^{^e) ■ du {^[g^ + \eg''^ (d.r) ® (9,r)^ | • Sr^/d^d^'^^ (Bl) 

= Jp' (Vd^ - 2g^^"(5,r) • (a,,5r) v/d^d"-| 

= - y {dj {p' (2Vd^ - l) .g'^(a,r) Vd^) } • -^r d"-| 

= - / { (V^'O ■ ^fe (p' (2v/d^ - l) 9''{drr) ® (5,r)) } • Sr v/d^d"^,^, (B2) 



which is summarized as 



with 



5 j (e-p' (v/d^ -l))dV^ j {(ye) ■ dk {P'H^^v) ® {d,v))} ■ 5v v/d^d"^| (B3) 

P'' = -S (g^ + ^£5'^) + P' (2 v/d^ti - 1) 5'^ = -^gj-' + (p' - E) g^^ (B4) 



where the last equal sign is due to v'detg = 1 . Then, rewriting the undetermined multiplier as p' = p + -E — 5, we 
obtain Eq. HI]). 
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According to Nakahara (private communication), a typ- 
ical value of the Reynolds number in such cases is Rh = 
100 on the basis of the layer thickness H, or Rl = 2000 
on the basis of the horizontal length scale L of the con- 
tainer. While Rh ~ 100 is typically not large enough 
to cause a transition to turbulence (in the usual sense 
of the word), we may expect a different kind of "turbu- 
lence" such as a two-dimensional chaotic flow maintained 
by horizontal forcing. 

The presence of the viscous drag is not explicitly stated in 
Duran's book [231, but it seems to be implicitly assumed 
by stating that x stops at the position satisfying kx — 
mG sin 6. 



